here::set_here()
## File .here already exists in /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/scripts
suppressPackageStartupMessages({
  library(slingshot)
  library(tradeSeq)
  library(SingleCellExperiment)
  library(cowplot)
  library(rgl)
  library(clusterExperiment)
  library(RColorBrewer)
  library(aggregation)
  library(ggplot2)
  library(pheatmap)
  library(wesanderson)
  library(UpSetR)
  library(gridExtra)
  library(NMF)
  library(nichenetr)
  library(Seurat)
  library(dplyr)
  library(tidyverse)
  library(circlize)
  library(ComplexHeatmap)
  library(NMF)
  library(msigdbr)
  library(fgsea)
})
colby <- function(values, g=4){
  if(is(values, "character")){
    cols <- as.numeric(as.factor(values))
    return(cols)
  }
  if(is(values, "factor")){
    ncl <- nlevels(values)
    if(ncl > 9){
          colpal <- c(RColorBrewer::brewer.pal(9, 'Set1'), wesanderson::wes_palette("Darjeeling1", n=ncl-9, type="continuous"))
    } else {
       colpal <- RColorBrewer::brewer.pal(9, 'Set1')
    }
    cols <- colpal[as.numeric(values)]
    return(cols)
  }
  if(is(values, "numeric")){
    pal <- wesanderson::wes_palette("Zissou1", n=g, type="continuous")
    gg <- Hmisc::cut2(values, g=g)
    if(nlevels(gg) < g){
      nl <- nlevels(gg)
      if(nl == 2) pal <- pal[c(1,g)]
    }
    cols <- pal[gg]
    return(cols)
  }
}

Import data

sds <- readRDS("../data/finalTrajectory/sling.rds")
counts <- readRDS("../data/finalTrajectory/counts_noResp_noMV.rds")
counts <- round(counts)
sce <- readRDS("../data/finalTrajectory/sce_tradeSeq20200904.rds")
load("../data/ALL_TF.Rda")
clDatta <- readRDS("../data/finalTrajectory/dattaCl_noResp_noMV.rds")

Trajectory

plot3d(reducedDim(sds), aspect = 'iso', col=brewer.pal(9,'Set1')[clDatta], alpha=.6)
plot3d.SlingshotDataSet(sds, add=TRUE, col=c("black"), lwd=4)

Nichenet analysis

library(nichenetr)
library(tidyverse)

Load ligand-target interaction matrix

ligand_target_matrix = readRDS(url("https://zenodo.org/record/3260758/files/ligand_target_matrix.rds"))
# convert human to mouse
colnames(ligand_target_matrix) = ligand_target_matrix %>% colnames() %>% convert_human_to_mouse_symbols()
rownames(ligand_target_matrix) = ligand_target_matrix %>% rownames() %>% convert_human_to_mouse_symbols()

Define expressed genes in sender and receiver population

Sender population is regenerated HBC. Receiver population is activated HBC. We define an expressed gene is a gene being expressed (non-zero) in at least 10% of the cells in that cell type.

expressed_genes_sender <- rownames(counts)[rowMeans(counts[, clDatta == "HBC"] > 0) > 0.05]
expressed_genes_receiver <- rownames(counts)[rowMeans(counts[, clDatta == "HBC*"] > 0) > 0.05]
length(expressed_genes_sender)
## [1] 8453
length(expressed_genes_receiver)
## [1] 8021

Nichenet analysis: HBC vs HBC* gene set of interest

Define gene set of interest through differential expression and background genes

clDatta <- droplevels(clDatta)
clDatta <- relevel(clDatta, ref="HBC*")
library(glmGamPoi)
fit <- glm_gp(counts, 
              design = model.matrix(~ clDatta))
deRes <- test_de(fit, contrast = "clDattaHBC")
deGenes <- deRes$name[deRes$adj_pval <= 0.01 & abs(deRes$lfc) > log2(2)]
geneset_oi <- deGenes

background_expressed_genes = expressed_genes_receiver %>% .[. %in% rownames(ligand_target_matrix)]

Define potential ligands

lr_network = readRDS(url("https://zenodo.org/record/3260758/files/lr_network.rds"))
lr_network = lr_network %>% mutate(from = convert_human_to_mouse_symbols(from), to = convert_human_to_mouse_symbols(to)) %>% drop_na()

ligands = lr_network %>% pull(from) %>% unique()
expressed_ligands = intersect(ligands,expressed_genes_sender)
receptors = lr_network %>% pull(to) %>% unique()
expressed_receptors = intersect(receptors,expressed_genes_receiver)

lr_network_expressed = lr_network %>% filter(from %in% expressed_ligands & to %in% expressed_receptors) 
head(lr_network_expressed)
## # A tibble: 6 x 4
##   from    to        source         database
##   <chr>   <chr>     <chr>          <chr>   
## 1 Il6     Il6st     kegg_cytokines kegg    
## 2 Ctf1    Lifr      kegg_cytokines kegg    
## 3 Ctf1    Il6st     kegg_cytokines kegg    
## 4 Tslp    Crlf2     kegg_cytokines kegg    
## 5 Tnfsf12 Tnfrsf12a kegg_cytokines kegg    
## 6 Tgfb1   Tgfbr1    kegg_cytokines kegg
potential_ligands = lr_network_expressed %>% pull(from) %>% unique()
head(potential_ligands)
## [1] "Il6"     "Ctf1"    "Tslp"    "Tnfsf12" "Tgfb1"   "Tgfb2"

Ligand activity analysis

ligand_activities = predict_ligand_activities(geneset = geneset_oi,
                                              background_expressed_genes = background_expressed_genes, 
                                              ligand_target_matrix = ligand_target_matrix, 
                                              potential_ligands = potential_ligands)

ligand_activities %>% arrange(-pearson)
## # A tibble: 132 x 4
##    test_ligand auroc  aupr  pearson
##    <chr>       <dbl> <dbl>    <dbl>
##  1 Hmgb2       0.506 0.326  0.151  
##  2 Dsc3        0.550 0.318  0.0824 
##  3 App         0.493 0.313  0.0718 
##  4 Apoe        0.478 0.270  0.0295 
##  5 Tgfb1       0.483 0.277  0.0288 
##  6 Tgfb2       0.469 0.275  0.0187 
##  7 Calr        0.499 0.280  0.0139 
##  8 Ltf         0.477 0.272  0.0114 
##  9 Il6         0.479 0.267  0.00132
## 10 Col4a1      0.505 0.274 -0.00148
## # … with 122 more rows
best_upstream_ligands = ligand_activities %>% top_n(20, pearson) %>% arrange(-pearson) %>% pull(test_ligand)

# Top 20 of most active ligands:
best_upstream_ligands
##  [1] "Hmgb2"  "Dsc3"   "App"    "Apoe"   "Tgfb1"  "Tgfb2"  "Calr"   "Ltf"   
##  [9] "Il6"    "Col4a1" "Jag1"   "Mif"    "Adam17" "Lamb2"  "Ctf1"   "Fgf1"  
## [17] "Sema5a" "Nrxn1"  "Has2"   "Lefty1"
# threshold
p_hist_lig_activity = ggplot(ligand_activities, aes(x=pearson)) + 
  geom_histogram(color="black", fill="darkorange")  + 
  # geom_density(alpha=.1, fill="orange") +
  geom_vline(aes(xintercept=min(ligand_activities %>% top_n(20, pearson) %>% pull(pearson))), color="red", linetype="dashed", size=1) + 
  labs(x="ligand activity (PCC)", y = "# ligands") +
  theme_classic()
p_hist_lig_activity
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Target genes of top ligands

active_ligand_target_links_df = best_upstream_ligands %>% lapply(get_weighted_ligand_target_links,geneset = geneset_oi, ligand_target_matrix = ligand_target_matrix, n = 250) %>% bind_rows()

active_ligand_target_links = prepare_ligand_target_visualization(ligand_target_df = active_ligand_target_links_df, ligand_target_matrix = ligand_target_matrix, cutoff = 0.25)
order_ligands = intersect(best_upstream_ligands, colnames(active_ligand_target_links)) %>% rev()
order_targets = intersect(active_ligand_target_links_df$target, rownames(active_ligand_target_links)) %>% unique()
vis_ligand_target = active_ligand_target_links[order_targets,order_ligands] %>% t()


p_ligand_target_network = vis_ligand_target %>% make_heatmap_ggplot("Prioritized HBC-ligands","HBC* target genes", color = "purple",legend_position = "top", x_axis_position = "top",legend_title = "Regulatory potential") + scale_fill_gradient2(low = "whitesmoke",  high = "purple", breaks = c(0,0.005,0.01)) + theme(axis.text.x = element_text(face = "italic"))
## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.
color <- colorRampPalette(rev(brewer.pal(n = 7, name = "RdYlBu")))(100)
breaks <- seq(0, 0.01, length=100)
breaks <- c(breaks, 0.05)
pheatmap::pheatmap(vis_ligand_target, cluster_rows=FALSE, cluster_cols=FALSE,
         col=color, breaks=breaks, border_color = NA)

Circos plot

current_plot_df <- vis_ligand_target
active_ligand_target_links_df = best_upstream_ligands %>% lapply(get_weighted_ligand_target_links,geneset = geneset_oi, ligand_target_matrix = ligand_target_matrix, n = 250) %>% bind_rows()

cutoff_include_all_ligands = active_ligand_target_links_df$weight %>% quantile(0.66)
active_ligand_target_links_df_circos = active_ligand_target_links_df %>% filter(weight > cutoff_include_all_ligands)
ligands_to_remove = setdiff(active_ligand_target_links_df$ligand %>% unique(), active_ligand_target_links_df_circos$ligand %>% unique())
targets_to_remove = setdiff(active_ligand_target_links_df$target %>% unique(), active_ligand_target_links_df_circos$target %>% unique())
circos_links = active_ligand_target_links_df %>% filter(!target %in% targets_to_remove &!ligand %in% ligands_to_remove)
logfc = deRes$lfc
logfc_filter = deRes$name[abs(logfc) >= 2]
circos_links = circos_links %>% filter(target %in% logfc_filter)
# logfc_sender = deRes$lfc[deRes$name %in% rownames(current_plot_df)]
# circos_links = circos_links %>% arrange(logfc[circos_links$target])
#     

cdm_res = chordDiagram(circos_links, annotationTrack = "grid", preAllocateTracks = 1)
## Note: The second link end is drawn out of sector 'Ccl12'.
## Note: The first link end is drawn out of sector 'Mif'.
abs_max = max(logfc, na.rm = TRUE)
    # https://jokergoo.github.io/circlize_book/book/legends.html#legends
    # https://jokergoo.github.io/circlize_book/book/a-complex-example-of-chord-diagram.html
col_fun = colorRamp2(c(-abs_max, 0, abs_max), c("dodgerblue4", "white", "red"))
col_fun2 = colorRamp2(c(0, max(ligand_activities[, "pearson"])), c("white", "green4"))
lgd_links = Legend(at = c(-4, -2, 0, 2, 4), col_fun = col_fun, title_position = "topleft", title = "LogFC")
lgd_act = Legend(at = c(0, ceiling(100*max(ligand_activities[, "pearson"])))/100, col_fun = col_fun2, title_position = "topleft", title = "Ligand Activity")
lgd_list_vertical = packLegend(lgd_links, lgd_act)

ylim = get.cell.meta.data("ylim", sector.index = circos_links$ligand[1], track.index = 1)
y1 = ylim[2] - (ylim[2] - ylim[1])*0.9
y2 = ylim[2] - (ylim[2] - ylim[1])*0.75
y3 = ylim[2] - (ylim[2] - ylim[1])*0.74
y4 = ylim[2] - (ylim[2] - ylim[1])*0.59
    
for(i in seq_len(nrow(cdm_res))) {
  if(cdm_res$value1[i] > 0) {
      circos.rect(
          cdm_res[i, "x1"], y1, 
          cdm_res[i, "x1"] - abs(cdm_res[i, "value1"]), y1 + (y2-y1)*0.7,
          col = col_fun(min(5, max(-5, logfc[cdm_res$rn[i]]))),
          border = col_fun(min(5, max(-5, logfc[cdm_res$rn[i]]))),
          sector.index = cdm_res$rn[i], track.index = 1)
      circos.rect(
          cdm_res[i, "x1"], y3 + (y4-y3)*0.3, 
          cdm_res[i, "x1"] - abs(cdm_res[i, "value1"]), y4,
          col = col_fun2(ligand_activities %>% filter(test_ligand == cdm_res$rn[i]) %>% .$pearson),
          border = col_fun2(ligand_activities %>% filter(test_ligand == cdm_res$rn[i]) %>% .$pearson),
          sector.index = cdm_res$rn[i], track.index = 1)
      circos.rect(
          cdm_res[i, "x2"], y1 + (y2-y1)*0.3, 
          cdm_res[i, "x2"] - abs(cdm_res[i, "value1"]), y2,
          col = col_fun(min(5, max(-5, logfc[cdm_res$cn[i]]))),
          border = col_fun(min(5, max(-5, logfc[cdm_res$cn[i]]))),
          sector.index = cdm_res$cn[i], track.index = 1)
        }
    }
## Note: 3 points are out of plotting region in sector 'Ccl12', track '1'.
## Note: 3 points are out of plotting region in sector 'Mif', track '1'.
## Note: 3 points are out of plotting region in sector 'Mif', track '1'.
    circos.trackPlotRegion(track.index = 1, panel.fun = function(x, y) {
      xlim = get.cell.meta.data("xlim")
      sector.name = get.cell.meta.data("sector.index")
      circos.text(mean(xlim), y4 + .1, sector.name, facing = "clockwise", niceFacing = TRUE, adj = c(0, 0.5), cex = 0.5)
    }, bg.border = NA)
    draw(lgd_list_vertical, x = unit(0.95, "npc"), y = unit(0.85, "npc"), just = c("top", "right"))

    circos.clear()

Nichenet analysis: HBC* vs all

Define gene set of interest through differential expression and background genes

clDatta <- droplevels(clDatta)
clDatta <- relevel(clDatta, ref="HBC*")
library(glmGamPoi)
fit <- glm_gp(counts, 
              design = model.matrix(~ clDatta))
L <- matrix(0, nrow=ncol(fit$Beta), ncol=1)
rownames(L) <- colnames(fit$Beta)
L[,1] <- 1/5
deRes <- test_de(fit, 
                 contrast = L)
deGenes <- deRes$name[deRes$adj_pval <= 0.01 & abs(deRes$lfc) > log2(2)]
length(deGenes)
## [1] 5623
geneset_oi <- deGenes

background_expressed_genes = expressed_genes_receiver %>% .[. %in% rownames(ligand_target_matrix)]

Define potential ligands

lr_network = readRDS(url("https://zenodo.org/record/3260758/files/lr_network.rds"))
lr_network = lr_network %>% mutate(from = convert_human_to_mouse_symbols(from), to = convert_human_to_mouse_symbols(to)) %>% drop_na()

ligands = lr_network %>% pull(from) %>% unique()
expressed_ligands = intersect(ligands,expressed_genes_sender)
receptors = lr_network %>% pull(to) %>% unique()
expressed_receptors = intersect(receptors,expressed_genes_receiver)

lr_network_expressed = lr_network %>% filter(from %in% expressed_ligands & to %in% expressed_receptors) 
head(lr_network_expressed)
## # A tibble: 6 x 4
##   from    to        source         database
##   <chr>   <chr>     <chr>          <chr>   
## 1 Il6     Il6st     kegg_cytokines kegg    
## 2 Ctf1    Lifr      kegg_cytokines kegg    
## 3 Ctf1    Il6st     kegg_cytokines kegg    
## 4 Tslp    Crlf2     kegg_cytokines kegg    
## 5 Tnfsf12 Tnfrsf12a kegg_cytokines kegg    
## 6 Tgfb1   Tgfbr1    kegg_cytokines kegg
potential_ligands = lr_network_expressed %>% pull(from) %>% unique()
head(potential_ligands)
## [1] "Il6"     "Ctf1"    "Tslp"    "Tnfsf12" "Tgfb1"   "Tgfb2"

Ligand activity analysis

ligand_activities = predict_ligand_activities(geneset = geneset_oi,
                                              background_expressed_genes = background_expressed_genes, 
                                              ligand_target_matrix = ligand_target_matrix, 
                                              potential_ligands = potential_ligands)

ligand_activities %>% arrange(-pearson)
## # A tibble: 132 x 4
##    test_ligand auroc  aupr pearson
##    <chr>       <dbl> <dbl>   <dbl>
##  1 Dsc3        0.579 0.575  0.148 
##  2 Hmgb2       0.531 0.535  0.0970
##  3 App         0.531 0.549  0.0930
##  4 Col4a1      0.546 0.530  0.0761
##  5 Apoe        0.520 0.519  0.0737
##  6 Lamb2       0.545 0.530  0.0731
##  7 Tgfb1       0.526 0.525  0.0721
##  8 Gas6        0.542 0.526  0.0678
##  9 Lefty1      0.534 0.518  0.0633
## 10 Adam17      0.532 0.520  0.0630
## # … with 122 more rows
best_upstream_ligands = ligand_activities %>% top_n(15, pearson) %>% arrange(-pearson) %>% pull(test_ligand)

# Top 20 of most active ligands:
best_upstream_ligands
##  [1] "Dsc3"   "Hmgb2"  "App"    "Col4a1" "Apoe"   "Lamb2"  "Tgfb1"  "Gas6"  
##  [9] "Lefty1" "Adam17" "Tslp"   "Sema5a" "Sema3a" "Has2"   "Celsr1"
# threshold
p_hist_lig_activity = ggplot(ligand_activities, aes(x=pearson)) + 
  geom_histogram(color="black", fill="darkorange")  + 
  # geom_density(alpha=.1, fill="orange") +
  geom_vline(aes(xintercept=min(ligand_activities %>% top_n(20, pearson) %>% pull(pearson))), color="red", linetype="dashed", size=1) + 
  labs(x="ligand activity (PCC)", y = "# ligands") +
  theme_classic()
p_hist_lig_activity
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Target genes of top ligands

active_ligand_target_links_df = best_upstream_ligands %>% lapply(get_weighted_ligand_target_links,geneset = geneset_oi, ligand_target_matrix = ligand_target_matrix, n = 250) %>% bind_rows()

active_ligand_target_links = prepare_ligand_target_visualization(ligand_target_df = active_ligand_target_links_df, ligand_target_matrix = ligand_target_matrix, cutoff = 0.25)
order_ligands = intersect(best_upstream_ligands, colnames(active_ligand_target_links)) %>% rev()
order_targets = intersect(active_ligand_target_links_df$target, rownames(active_ligand_target_links)) %>% unique()
vis_ligand_target = active_ligand_target_links[order_targets,order_ligands] %>% t()


p_ligand_target_network = vis_ligand_target %>% make_heatmap_ggplot("Prioritized ligands","HBC* target genes", color = "purple",legend_position = "top", x_axis_position = "top",legend_title = "Regulatory potential") + scale_fill_gradient2(low = "whitesmoke",  high = "purple", breaks = c(0,0.005,0.01)) + theme(axis.text.x = element_text(face = "italic"))
## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.
color <- colorRampPalette(rev(brewer.pal(n = 7, name = "RdYlBu")))(100)
breaks <- seq(0, 0.01, length=100)
breaks <- c(breaks, 0.05)
pheatmap::pheatmap(vis_ligand_target, cluster_rows=FALSE, cluster_cols=FALSE,
         col=color, breaks=breaks, border_color = NA)

Circos plot

current_plot_df <- vis_ligand_target
active_ligand_target_links_df = best_upstream_ligands %>% lapply(get_weighted_ligand_target_links,geneset = geneset_oi, ligand_target_matrix = ligand_target_matrix, n = 250) %>% bind_rows()

cutoff_include_all_ligands = active_ligand_target_links_df$weight %>% quantile(0.66)
active_ligand_target_links_df_circos = active_ligand_target_links_df %>% filter(weight > cutoff_include_all_ligands)
ligands_to_remove = setdiff(active_ligand_target_links_df$ligand %>% unique(), active_ligand_target_links_df_circos$ligand %>% unique())
targets_to_remove = setdiff(active_ligand_target_links_df$target %>% unique(), active_ligand_target_links_df_circos$target %>% unique())
circos_links = active_ligand_target_links_df %>% filter(!target %in% targets_to_remove &!ligand %in% ligands_to_remove)
logfc = deRes$lfc
logfc_filter = deRes$name[abs(logfc) >= 2]
circos_links = circos_links %>% filter(target %in% logfc_filter)
# logfc_sender = deRes$lfc[deRes$name %in% rownames(current_plot_df)]
# circos_links = circos_links %>% arrange(logfc[circos_links$target])
#     

cdm_res = chordDiagram(circos_links, annotationTrack = "grid", preAllocateTracks = 1)
## Note: The first link end is drawn out of sector 'Hmgb2'.
## Note: The second link end is drawn out of sector 'Ets1'.
## Note: The second link end is drawn out of sector 'Cd86'.
## Note: The first link end is drawn out of sector 'Tslp'.
abs_max = max(logfc, na.rm = TRUE)
    # https://jokergoo.github.io/circlize_book/book/legends.html#legends
    # https://jokergoo.github.io/circlize_book/book/a-complex-example-of-chord-diagram.html
col_fun = colorRamp2(c(-abs_max, 0, abs_max), c("dodgerblue4", "white", "red"))
col_fun2 = colorRamp2(c(0, max(ligand_activities[, "pearson"])), c("white", "green4"))
lgd_links = Legend(at = c(-4, -2, 0, 2, 4), col_fun = col_fun, title_position = "topleft", title = "LogFC")
lgd_act = Legend(at = c(0, ceiling(100*max(ligand_activities[, "pearson"])))/100, col_fun = col_fun2, title_position = "topleft", title = "Ligand Activity")
lgd_list_vertical = packLegend(lgd_links, lgd_act)

ylim = get.cell.meta.data("ylim", sector.index = circos_links$ligand[1], track.index = 1)
y1 = ylim[2] - (ylim[2] - ylim[1])*0.9
y2 = ylim[2] - (ylim[2] - ylim[1])*0.75
y3 = ylim[2] - (ylim[2] - ylim[1])*0.74
y4 = ylim[2] - (ylim[2] - ylim[1])*0.59
    
for(i in seq_len(nrow(cdm_res))) {
  if(cdm_res$value1[i] > 0) {
      circos.rect(
          cdm_res[i, "x1"], y1, 
          cdm_res[i, "x1"] - abs(cdm_res[i, "value1"]), y1 + (y2-y1)*0.7,
          col = col_fun(min(5, max(-5, logfc[cdm_res$rn[i]]))),
          border = col_fun(min(5, max(-5, logfc[cdm_res$rn[i]]))),
          sector.index = cdm_res$rn[i], track.index = 1)
      circos.rect(
          cdm_res[i, "x1"], y3 + (y4-y3)*0.3, 
          cdm_res[i, "x1"] - abs(cdm_res[i, "value1"]), y4,
          col = col_fun2(ligand_activities %>% filter(test_ligand == cdm_res$rn[i]) %>% .$pearson),
          border = col_fun2(ligand_activities %>% filter(test_ligand == cdm_res$rn[i]) %>% .$pearson),
          sector.index = cdm_res$rn[i], track.index = 1)
      circos.rect(
          cdm_res[i, "x2"], y1 + (y2-y1)*0.3, 
          cdm_res[i, "x2"] - abs(cdm_res[i, "value1"]), y2,
          col = col_fun(min(5, max(-5, logfc[cdm_res$cn[i]]))),
          border = col_fun(min(5, max(-5, logfc[cdm_res$cn[i]]))),
          sector.index = cdm_res$cn[i], track.index = 1)
        }
    }
## Note: 3 points are out of plotting region in sector 'Hmgb2', track '1'.
## Note: 3 points are out of plotting region in sector 'Hmgb2', track '1'.
## Note: 3 points are out of plotting region in sector 'Ets1', track '1'.
## Note: 3 points are out of plotting region in sector 'Cd86', track '1'.
## Note: 3 points are out of plotting region in sector 'Tslp', track '1'.
## Note: 3 points are out of plotting region in sector 'Tslp', track '1'.
    circos.trackPlotRegion(track.index = 1, panel.fun = function(x, y) {
      xlim = get.cell.meta.data("xlim")
      sector.name = get.cell.meta.data("sector.index")
      circos.text(mean(xlim), y4 + .1, sector.name, facing = "clockwise", niceFacing = TRUE, adj = c(0, 0.5), cex = 0.5)
    }, bg.border = NA)
    draw(lgd_list_vertical, x = unit(0.95, "npc"), y = unit(0.85, "npc"), just = c("top", "right"))

    circos.clear()

Random input in terms of DE genes

prepare_ligand_target_visualization_kvdb <- function(ligand_target_df,
                                                     ligand_target_matrix, 
                                                     cutoff = 0.25) 
{
    cutoff_include_all_ligands = ligand_target_df$weight %>% 
        quantile(cutoff, na.rm=TRUE)
    ligand_target_matrix_oi = ligand_target_matrix
    ligand_target_matrix_oi[ligand_target_matrix_oi < cutoff_include_all_ligands] = 0
    
    ligand_target_vis = ligand_target_matrix_oi[ligand_target_df$target %>% 
        unique(), ligand_target_df$ligand %>% unique()]
    dim(ligand_target_vis) = c(length(ligand_target_df$target %>% 
        unique()), length(ligand_target_df$ligand %>% unique()))
    all_targets = ligand_target_df$target %>% unique()
    all_ligands = ligand_target_df$ligand %>% unique()
    rownames(ligand_target_vis) = all_targets
    colnames(ligand_target_vis) = all_ligands
    keep_targets = all_targets[ligand_target_vis %>% apply(1, 
        sum) > 0]
    keep_ligands = all_ligands[ligand_target_vis %>% apply(2, 
        sum) > 0]
    ligand_target_vis_filtered = ligand_target_vis[keep_targets, 
        keep_ligands]
    if (is.matrix(ligand_target_vis_filtered)) {
        rownames(ligand_target_vis_filtered) = keep_targets
        colnames(ligand_target_vis_filtered) = keep_ligands
    }
    else {
        dim(ligand_target_vis_filtered) = c(length(keep_targets), 
            length(keep_ligands))
        rownames(ligand_target_vis_filtered) = keep_targets
        colnames(ligand_target_vis_filtered) = keep_ligands
    }
    if (nrow(ligand_target_vis_filtered) > 1 & ncol(ligand_target_vis_filtered) > 
        1) {
        distoi = dist(1 - cor(t(ligand_target_vis_filtered)))
        hclust_obj = hclust(distoi, method = "ward.D2")
        order_targets = hclust_obj$labels[hclust_obj$order]
        distoi_targets = dist(1 - cor(ligand_target_vis_filtered))
        hclust_obj = hclust(distoi_targets, method = "ward.D2")
        order_ligands = hclust_obj$labels[hclust_obj$order]
    }
    else {
        order_targets = rownames(ligand_target_vis_filtered)
        order_ligands = colnames(ligand_target_vis_filtered)
    }
    vis_ligand_target_network = ligand_target_vis_filtered[order_targets, 
        order_ligands]
    dim(vis_ligand_target_network) = c(length(order_targets), 
        length(order_ligands))
    rownames(vis_ligand_target_network) = order_targets
    colnames(vis_ligand_target_network) = order_ligands
    return(vis_ligand_target_network)
}

runNicheNet <- function(genesOfInterest, nTopLigands,
                        expressed_genes_receiver,
                        expressed_genes_sender){
  
  
  geneset_oi <- genesOfInterest
  background_expressed_genes = expressed_genes_receiver %>% .[. %in% rownames(ligand_target_matrix)]
  
  lr_network = readRDS(url("https://zenodo.org/record/3260758/files/lr_network.rds"))
  lr_network = lr_network %>% mutate(from = convert_human_to_mouse_symbols(from), to = convert_human_to_mouse_symbols(to)) %>% drop_na()
  
  ligands = lr_network %>% pull(from) %>% unique()
  expressed_ligands = intersect(ligands,expressed_genes_sender)
  receptors = lr_network %>% pull(to) %>% unique()
  expressed_receptors = intersect(receptors,expressed_genes_receiver)
  
  lr_network_expressed = lr_network %>% filter(from %in% expressed_ligands & to %in% expressed_receptors) 
  head(lr_network_expressed)
  
  potential_ligands = lr_network_expressed %>% pull(from) %>% unique()
  head(potential_ligands)
  
  ligand_activities = predict_ligand_activities(geneset = geneset_oi,
                                              background_expressed_genes = background_expressed_genes, 
                                              ligand_target_matrix = ligand_target_matrix, 
                                              potential_ligands = potential_ligands)

  ligand_activities %>% arrange(-pearson)
  best_upstream_ligands = ligand_activities %>% top_n(nTopLigands, pearson) %>% arrange(-pearson) %>% pull(test_ligand)
  
  # Top 20 of most active ligands:
  best_upstream_ligands
  
  # threshold
  p_hist_lig_activity = ggplot(ligand_activities, aes(x=pearson)) + 
    geom_histogram(color="black", fill="darkorange")  + 
    # geom_density(alpha=.1, fill="orange") +
    geom_vline(aes(xintercept=min(ligand_activities %>% top_n(nTopLigands, pearson) %>% pull(pearson))), color="red", linetype="dashed", size=1) + 
    labs(x="ligand activity (PCC)", y = "# ligands") +
    theme_classic()
  p_hist_lig_activity
  
  ### changed
      # active_ligand_target_links_df = best_upstream_ligands %>% lapply(get_weighted_ligand_target_links,geneset = geneset_oi, ligand_target_matrix = ligand_target_matrix, n = 250) %>% bind_rows()

  
    active_ligand_target_links_df = best_upstream_ligands %>%
      lapply(get_weighted_ligand_target_links,
             geneset = geneset_oi, 
             ligand_target_matrix = ligand_target_matrix, 
             n = 250) %>% 
      bind_rows()
    # added:
    active_ligand_target_links_df <- active_ligand_target_links_df[!is.na(active_ligand_target_links_df$target),]
  
  active_ligand_target_links <- prepare_ligand_target_visualization_kvdb(
    ligand_target_df = active_ligand_target_links_df, 
    ligand_target_matrix = ligand_target_matrix, 
    cutoff = 0.25)
  order_ligands = intersect(best_upstream_ligands, colnames(active_ligand_target_links)) %>% rev()
  order_targets = intersect(active_ligand_target_links_df$target, rownames(active_ligand_target_links)) %>% unique()
  vis_ligand_target = active_ligand_target_links[order_targets,order_ligands] %>% t()
  
  
  p_ligand_target_network = vis_ligand_target %>% make_heatmap_ggplot("Prioritized ligands","HBC* target genes", color = "purple",legend_position = "top", x_axis_position = "top",legend_title = "Regulatory potential") + scale_fill_gradient2(low = "whitesmoke",  high = "purple", breaks = c(0,0.005,0.01)) + theme(axis.text.x = element_text(face = "italic"))
  
  
  color <- colorRampPalette(rev(brewer.pal(n = 7, name = "RdYlBu")))(100)
  breaks <- seq(0, 0.01, length=100)
  breaks <- c(breaks, 0.05)
  ph <- pheatmap::pheatmap(vis_ligand_target, cluster_rows=FALSE, cluster_cols=FALSE,
           col=color, breaks=breaks, border_color = NA)
  return(ph)
}

5000 genes as input

length(deGenes) ## 5623
## [1] 5623
phList <- list()
for(kk in 1:5){
  set.seed(kk)
  genes <- sample(rownames(counts), 5000)
  phList[[kk]] <- runNicheNet(genesOfInterest=genes, nTopLigands=20,
                        expressed_genes_receiver,
                        expressed_genes_sender)
  print(phList[[kk]])
}
## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.
## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.

## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.

## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.

## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.

2000 genes as input

phList2k <- list()
for(kk in 1:3){
  set.seed(kk)
  genes <- sample(rownames(counts), 2000)
  phList2k[[kk]] <- runNicheNet(genesOfInterest=genes, nTopLigands=20,
                        expressed_genes_receiver,
                        expressed_genes_sender)
  print(phList2k[[kk]])
}
## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.
## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.

## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.

200 genes as input

phList2k <- list()
for(kk in 1:6){
  set.seed(kk)
  genes <- sample(rownames(counts), 200)
  phList2k[[kk]] <- runNicheNet(genesOfInterest=genes, nTopLigands=15,
                        expressed_genes_receiver,
                        expressed_genes_sender)
  print(phList2k[[kk]])
}
## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.
## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.

## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.

## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.

## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.

## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.